arXiv:1508.04788vl [cond-mat.mtrl-sci] 19 Aug 2015 


The Nature of the Interlayer Interaction in Bulk and Few-Layer Phosphorus 


L. Shulenburger , 1 A.D. Baczewski , 1 Z. Zhu , 2 J. Guan , 2 and D. Tomanek 2 

1 Sandia National Laboratories, Albuquerque NM 87185 
2 Michigan State University, East Lansing MI 48825 

An outstanding challenge of theoretical electronic structure is the description of van der Waals 
(vdW) interactions in molecules and solids. Renewed interest in resolving this is in part motivated 
by the technological promise of layered systems including graphite, transition metal dichalcogenides, 
and more recently, black phosphorus, in which the interlayer interaction is widely believed to be 
dominated by these types of forces. We report a series of quantum Monte Carlo (QMC) calculations 
for bulk black phosphorus and related few-layer phosphorene, which elucidate the nature of the 
forces that bind these systems and provide benchmark data for the energetics of these systems. We 
find a significant charge redistribution due to the interaction between electrons on adjacent layers. 
Comparison to density functional theory (DFT) calculations indicate not only wide variability even 
among different vdW corrected functionals, but the failure of these functionals to capture the trend 
of reorganization predicted by QMC. The delicate interplay of steric and dispersive forces between 
layers indicate that few-layer phosphorene presents an unexpected challenge for the development of 
vdW corrected DFT. 


There is growing interest in understanding and cor¬ 
rectly describing the nature of the weak interlayer in¬ 
teraction in layered systems ranging from few-layer 
graphene 1 to transition metal dichalcogenidePI such as 
M0S2 and few-layer phosphorend^, which display unique 
electronic properties and bear promise for device applica¬ 
tions. Anticipating that challenges concerning the stabil¬ 
ity and isolation of single- to few-layer phosphorene can 
be overcomd^, this system with its unique electronic 3 4 
and opticaP properties is attracting particular interest. 
It displays a high and anisotropic carrier mobility 7 8 and 
a robust band gap that depends sensitively on the in-layer 
strain 3 . Progress in device fabrication? 9 indicates clearly 
that phosphorene holds technological promise. Since the 
fundamental band gap depends sensitively on the num¬ 
ber of layer J 3 , understanding the nature of the interlayer 
interaction is particulary important. 

The standard approach to describe the interlayer inter¬ 
action in layered solids has been based on DFT. Whereas 
DFT is - in principle - capable of describing the total en¬ 
ergy of any system in the ground state exactly, current 
implementations describe the effects of electron exchange 
and correlation only in an approximate manner. In most 
covalent and ionic solids of interest, the specific treat¬ 
ment of exchange and correlation of electrons does not 
play a crucial role and commonly used local or semi-local 
exchange-correlation functionals of the electron density 
are adequate. This approach may, however, not be ade¬ 
quate in complex and weakly bonded systems 10 including 
black phosphorus 11 . This is illustrated in Fig. [TJ which 
displays large differences between interlayer interaction 
energies in bulk and bilayer black phosphorus, obtained 
using different DFT functionals contained in the VASP 
software packagd^ 14 . The large spread of the interlayer 
energies predicted using these functionals, some of which 
include van der Waals (vdW) corrections, illustrates the 
gravity of the issue. 

Even though vdW corrected exchange-correlation 
functionals can improve the predicted geometry of lay- 
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FIG. 1. Binding energy per atom A E as a function of the 
interlayer spacing d in AB stacked (A) bulk and (B) bi¬ 
layer black phosphorus. QMC results obtained using dif¬ 
fusion Monte Carlo (DMC) are compared to DFT with 
LDAP, PBEP, AMQJpl DFT+DlP 1 , vdW-DFiP 1 , vdW- 
optB86b 20 21 and HSEOfP^ exchange-correlation functionals. 
The lines connecting the data points are Morse fits that ex¬ 
trapolate to A E — 0 for d— 00 . The vertical dashed line in¬ 
dicates the observed interlayer spacing d e (expt.) = 5.2365 A 
in the bulk structure in (A) and the optimum value based on 
DMC in the bilayer in (B). Side views of the geometries are 
shown in the insets. 


ered phosphorus, other quantities such as the transition 
pressure from bulk orthorhombic to rhombohedral phases 
are significantly underestimated 1 ^. This comes as no sur¬ 
prise, since recent benchmarks in dense hydrogen indi¬ 
cate that variations in the relative accuracy are observed 
among numerous vdW corrected functionals depending 
on the quantity being calculated 23 . Even though the na¬ 
ture of the interactions in phosphorus is quite different 
from hydrogen, the difficulty to properly balance vdW 
and other interactions likely persists. 

A superior way to obtain insight into the nature of the 
interlayer interaction requires a computational approach 
that treats electron exchange and correlation adequately 
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and on the same footing as covalent and ionic interac¬ 
tions. Unlike in DFT, where the fundamental quantity is 
the electron density, the fundamental quantity in QMC 
calculations is the properly antisymmetrized all-valence- 
electron wavefunction that explicitly describes the cor¬ 
relation of electrons. Therefore, the weak interlayer in¬ 
teraction in layered systems obtained using QMC is ex¬ 
pected to be more precise than the DFT-based counter¬ 
part, and the electron density obtained by QMC likely 
provides a better representation of the true charge den¬ 
sity than DFT. Should the QMC-based electron distri¬ 
bution in a few-layer system differ in a significant and 
nontrivial manner from a superposition of electron den¬ 
sities in isolated monolayers, we may surely conclude that 
the interlayer interaction in such a system is not purely 
dispersive in nature. 

In the last decade, QMC methods have demonstrated 
considerable promise as a high accuracy first principles 
method for studying solids 24 and are especially well- 
suited to vdW bound systems 25 31 . On one hand, the 
absence of approximations in treating the interaction be¬ 
tween electrons makes it an ideal method for layered ma¬ 
terials with a competition between different types of in¬ 
teractions. On the other hand, the computational cost 
associated QMC calculations is typically 1,000 — 10,000 
times higher than that of comparable semilocal DFT cal¬ 
culations. Even though this cost is mitigated by its su¬ 
perior parallel scalability, QMC is - for the time being - 
most useful for benchmark calculations. 

We have performed QMC calculations using the diffu¬ 
sion Monte Carlo (DMC) approach, as described in the 
Supplemental Material, in which details of the pseudopo¬ 
tential and an intensive procedure for converging finite 
size effects are included. One of the key approximations 
in our calculation is bias due to a fixed nodal surface. 
This was not anticipated to be significant in black phos¬ 
phorus as the binding of interest occurs in a region of low 
electronic density for which the degree of nodal nonlin¬ 
earity is expected to be low 32 i Nevertheless, some DMC 
calculations were carried out using orbitals from LDA, 
PBE, and vdW-optb86B functionals, to investigate the 
impact of the fixed node approximation. As the nodal 
surface associated with LDA orbitals were found to give 
the lowest energy and our method is variational, the LDA 
orbitals were subsequently used in all cases. 

Our DMC results for the interlayer binding curves in 
bulk and bilayer black phosphorus are shown by the solid 
lines in Fig. [lj The total energy difference between the 
monolayer and the bulk system indicates that the bind¬ 
ing energy of phosphorene sheets in black phosphorus 
is 81±6 meV/atom, which translates to a cleavage en¬ 
ergy of 22.4T1.6 meV per A 2 of the interface area. This 
is larger than that for many other layered materials 33 , 
but weak enough to allow mechanical exfoliation that 
has been reported. 3 4 

Comparing the DMC results to a variety of different 
DFT functionals, significant variability is evident. Re¬ 
sults are in agreement with intuition for LDA, PBE, and 


HSE06 which do not treat vdW explicitly - namely that 
LDA overbinds and reduces the interlayer spacing in com¬ 
parison to experiment, whereas GGA and HSE06 under¬ 
bind. The gradient corrected AM05, designed to prop¬ 
erly capture the energetics of the Airy gas at a jellium 
surface, has none of the spurious self-interaction in low 
density regions that cause LDA to overbind materials 
with vdW interactions. Nevertheless, it reproduces the 
bulk interlayer spacing correctly, but with only a quarter 
of the interlayer interaction energy predicted by DMC. 
That a vdW-free functional binds at all gives some indi¬ 
cation that the character of the interlayer interaction in 
these systems is not strictly vdW. 

From the wide array of available vdW functionals, 
we have selected 3 exemplary of increasing levels of 
sophistication. DFT+D2 is based upon an empirical 
correction to the total energy in the form of a simple 
pair-wise interaction parameterized by atomic Cq coeffi¬ 
cients. vdW-optB86b and vdW-DF2 are more advanced, 
both based upon improvements to the non-local vdW- 
DF functional 34 . Significant variability is evident in both 
energetics and the equilibrium interlayer spacing among 
these functionals. It is interesting to note that the least 
sophisticated of these functionals (DFT+D2) performs 
the best relative to DMC, as well as experiment in the 
case of the bulk system. 

Contrasting the bulk and bilayer binding curves, it is 
evident that DMC predicts that the interlayer interac¬ 
tion is not strictly additive. For an additive interaction, 
we should expect that the binding energy of the bilayer 
would be approximately 1/2 that of the bulk system be¬ 
cause both layers are missing half of their neighboring 
layers. Instead, DMC predicts that this is not the case 
and that the bilayer binding energy is instead 3/8 that 
of the bulk system. The results for vdW corrected DFT 
(DFT+D2, vdW-opB86b, and vdW-DF2) are more in¬ 
dicative of an additive interlayer interaction, giving us 
another indication that the nature of the interlayer bind¬ 
ing in black phosphorus is richer than a simple vdW in¬ 
teraction. 

A more direct indication of the character of the in¬ 
terlayer binding is the charge density difference induced 
by assembling the bulk system from isolated monolay¬ 
ers. We computed the quantity A p = p tot (bulk ) — 
22 ptot( mon °l a y ers ) using both DMC and DFT to in¬ 
vestigate this. The /i-norm of Ap over the unit cell is an 
indicator of the number of electrons being redistributed 
due to interlayer interaction. In all DFT functionals con¬ 
sidered in this study, this metric indicates a motion of 
fewer than 0.03 electrons per atom. In contrast, it pre¬ 
dicts a motion of 0.15 electrons per atom in DMC. To 
provide insight into the nature of this significant redis¬ 
tribution, the density difference is visualized in Fig. [2]4. 
Inspection indicates that charge is pushed out of the re¬ 
gion between layers and into the covalent bonds within 
each layer. This picture is well supported by basic chem¬ 
ical intuition. In an isolated layer, each atom is 3-fold 
coordinated with sp 3 bonding character and a single lone 
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pair protruding away from the layer. Bringing layers to¬ 
gether will increase the overlap between these lone pairs 
on adjacent layers and steric forces will tend to drive the 
affiliated lone pair charge closer to the layer on which it 
originated. 



FIG. 2. Electron density difference A p m p to t(bulk ) — 
^2 Ptot (monolayers) representing the charge redistribution 
caused by assembling the bulk structure from isolated mono- 
layers. (A) DMC isosurfaces bounding regions of excess 
electron density (dark brown) and electron deficiency (light 
brown), with respective values ±6.5 x 10 -3 e/A 3 . (B) (A p(z)) 
for DMC and several DFT functionals averaged across the 
x — y plane of the layers, with zfc indicating the relative po¬ 
sition of the plane in the unit cell. 

To further elucidate the role that this charge redistri¬ 
bution plays in the interlayer binding, the planar average 
of the charge density difference along planes perpendicu¬ 
lar to the interlayer axis for both DMC and several DFT 
functionals are illustrated in Fig. [2^. It remains evident 
that DMC predicts an average depletion of charge be¬ 
tween layers and commensurate accumulation in-layer. 
However, with the exception of vdW-DF2, this trend is 
not evident in any of the DFT calculations. Instead, DFT 
predicts a weak accumulation of charge between layers, 
indicating a minor covalent character in addition to what¬ 
ever vdW binding is present. This likely explains why a 
vdW-free functional like AM05 weakly binds the system 
near the correct interlayer spacing. The case of vdW- 
DF2 is of particular interest because it is the only func¬ 
tional that predicts a qualitatively similar, but consid¬ 
erably weaker, trend as DMC. Among DFT±D2, vdW- 
optB86b, and vdW-DF2, vdW-DF2 is the most theoreti¬ 
cally sophisticated so that it gets closer to the trend pre¬ 
dicted by DMC is unsurprising. Even so, that DFT±D2 
gives the best performance in terms of energetics and ge¬ 
ometry indicates that this functional may be getting the 
right answer for reasons that are not necessarily consis¬ 
tent with the many-body physics more explicitly explored 
through DMC. 

Recent work on self-consistent vdW functionals indi¬ 
cate that the charge redistribution induced by vdW in¬ 
teractions can play an important role in the energetics of 
highly polarizable systems 35 . In this work, the authors 


note that this subtle physics is consistent with an early 
observation by Feynman 36 in which the vdW interaction 
can be viewed as arising from an attractive interaction 
induced through a small accumulation of charge density 
between two mutually perturbed neutral systems. In the 
case of black phosphorus, we find that this picture of the 
vdW interaction is balanced by the steric redistribution 
of charge away from the region between layers. Based 
upon the results elucidated by our DMC calculations, we 
anticipate that getting this balance right may be a criti¬ 
cal requirement to examine in developing more advanced 
DFT functionals for layered compounds. 

The corresponding difference between the charge den¬ 
sity in an isolated monolayer and within a layer of bulk 
phosphorus is likely to affect the in-plane bonding and 
geometry. To see if this is indeed the case, we have cal¬ 
culated the energy change A E as a function of a stretch 
applied along the softer axis a\ of the sheets. Our results 
for the bulk system, presented in Fig. [3j indicate an ex¬ 
cellent agreement (to within half a percent) between the 
optimized lattice constant ai(theory ) = 4.404±0.019 A 
and the observed value 371 ai(expt.) = 4.374 A in the bulk 
structure. Further, we can see precisely how soft this axis 
is in both systems, with a deformation by |Aai|<0.3 A 
requiring an energy investment of only ^5 meV/atom in 
a monolayer or in bulk black phosphorus. This energy 
corresponds to a thermal energy of 60 K, and we should 
expect significant thermal fluctuations of the geometry of 
unsupported phosphorene sheets at ambient temperature 
and pressure. Most important, however, is the compari¬ 
son between a\ in the isolated monolayer and in the bulk 
structure. Our DMC results indicate a change in the in¬ 
plane stiffness along the soft axis and an ^2% reduction 
in the equilibrium lattice constant a\ in a monolayer from 
the bulk value. This is another indication of a charge re¬ 
distribution during the formation of a layered bulk struc¬ 
ture from monolayers that modifies the covalent interac¬ 
tion within the layers. This again supports our finding 
that the interlayer interaction in black phosphorus is not 
purely dispersive. 

To gain additional insight into the nature of the in¬ 
terlayer interaction, we compare in Fig. [4] the bonding 
within an AA and AB stacked bilayer as a function of 
the interlayer separation d. Our DMC results in Fig. [4] 
indicate that the AB stacking, which occurs in the bulk 
material, persists also in the bilayer. 

The cleavage energy of an AB stacked bilayer is 
16.6±2.2 meV/A 2 , thus 26% smaller than the bulk cleav¬ 
age energy of 22.4±1.6 meV/A 2 . The exfoliation energy 
associated with removing the topmost layer from the sur¬ 
face is expected to lie between these two values. The 
findings appear plausible also in view of the fact that 
in graphite, the cleavage energy is estimated to be 18% 
larger than the exfoliation energ}®. Though we expect 
the interlayer interaction to be mediated primarily by 
the 7r electrons in graphite and sp 3 -like lone pairs with 
a different character in black phosphorus, the ratio of 
the exfoliation and cleavage energy in the two systems 
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FIG. 3. DMC results for the relative total energy per atom 
A E as a function of the in-layer lattice constant a\ in the 
soft direction of a phosphorene monolayer and of bulk black 
phosphorus. The lines connecting the data points are Morse 
fits that extrapolate to A E — 0 for d—>• oo. The optimum 
lattice constant values for both structures are indicated by 
the vertical dashed lines. The observed value ai, e (expt.) = 
4.374 A in the bulk structure is indicated by the arrow. The 
monolayer geometry is shown in the inset. 


appears to be similar. 

The interlayer spacing d e = 5.272T0.023 A in the 
AB-stacked bilayer is about 1% larger than our calcu¬ 
lated value for the bulk material. In the less favorable 
AA stacking geometry, the binding energy is three times 
smaller than in the AB geometry, which should effec¬ 
tively prevent formation of stacking faults, at least in 
the absence of impurities. The significant difference be¬ 
tween the interaction in the AA and AB stacked bilayer 
suggests that the interaction between the sheets is more 
complicated than the vdW interaction between two ho¬ 
mogeneous slabs. Were the interlayer interaction purely 
dispersive, the registry of layers would not matter much 
and the AA and AB binding energies as well as interlayer 
separations should be nearly identical. 

To shed some light on the sensitivity of the interlayer 
bonding on the stacking sequence, we investigated the 
change in the charge density A p induced by the inter¬ 
action. Our results for AB and AA stackigns are shown 
in Fig. [4] C. Similar to the corresponding results for bulk 
black phosphorus in Fig. |2]B, the A p plots for the bilayer 
show a significant rearrangement of the electronic charge, 
a total of ^0.075 electrons per atom in both cases. In 
the AB-stacked bilayer, similar to the bulk system, we 
observe a depletion of the electron density in the region 
between the sheets and electron accumulation within the 
layers. The charge redistribution in the AA bilayer is 
significantly different, even including small regions in the 
interlayer space where the charge density increases. The 
large difference between Ap in the A A and AB stacked 
bilayer is inconsistent with purely dispersive bonding and 
explains why the the interlayer spacing and binding en¬ 
ergy are so different in the two systems. 



FIG. 4. DMC results for bonding in A A and AB stacked phos¬ 
phorene bilayers. (A) Geometry of an AA and AB stacked 
bilayer in top and side view. (B) DMC results for the rela¬ 
tive total energy per atom A E as a function of the interlayer 
spacing d. The lines connecting the data points are Morse fits 
that extrapolate to A E = 0 for d —> oo. (C) Electron den¬ 
sity differences for the AB and AA bilayers illustrated using 
isosurfaces and planar averaging as in Fig. [ 3 ] (with the same 
color coding). 


In summary, we studied the nature of the interlayer 
interaction in layered black phosphorus using quantum 
Monte Carlo calculations, which describe the correlation 
of electrons explicitly by an antisymmetrized all-valence- 
electron wavefunction and treat covalent and dispersive 
interactions on the same footing. Unlike in true vdW 
systems, we find that the interlayer interaction in few- 
layer phosphorene is associated with a significant charge 
redistribution between the in-layer and interlayer region, 
caused by changes in the non-local correlation of elec¬ 
trons in adjacent layers. Consequently, the resulting in¬ 
terlayer interaction can not be described properly by den¬ 
sity functional theory (DFT) augmented by mere semi¬ 
local vdW correction terms, and thus the designation 
‘van der Waals solids’ is improper for systems includ¬ 
ing few-layer phosphorene. Our results may be used as 
benchmarks for developing more sophisticated DFT func¬ 
tionals that should provide an improved description of 
non-local electron correlations in layered systems. 
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I. CONVERGENCE OF FINITE-SIZE EFFECTS 

One- and two-body finite size effects must be considered in both the bulk and planar periodic structures. One- 
body effects are corrected using canonical twist averaging®', while two-body effects are extrapolated. Controlling the 
two-body effects present a challenge for extrapolation consistent with reports of calculations done on graphite. 87 

For the bulk systems, we consider two types of tilings for generating these supercells. We expect finite size effects 
to converge more quickly when increasing the effective system size along the layering axis rather than in-plane. 
Consequently, we study two-body effects for 2x2x1, 3x3x1, 4x4x1, and 5x5x1 tilings (single layer) and 2x2x2, 3x3x2, 
and 4xe4x2 (bilayer) in the bulk. 

However, the peculiar nature of the electron correlation in two dimensional systems such as black phosphorus coupled 
with the meV accuracy required for this study requires a slightly different procedure for assessing and reducing the 
two body finite size effects introduced by simulating a supercell with periodic boundary conditions. The typically 
used procedures for this involve either utilizing a model interaction that removes spurious electron correlation, 88 
analyzing the behavior of the structure factor and two body jastrow factor for small values or momentum transfer 89 
or utilizin g ca lculations with density functionals designed to mimic the energetics of the electron correlation in finite 
supercells! 810 As generally practiced, all of these procedures rely on the assumption that the electron correlation is 
isotropic, a condition that is grossly violated in a layered material like black phosphorus. 



cr 1 (i/A) 


FIG. SI. Energy per atom of black phosphorus as a function of the inverse of the shortest distance between periodic images 
of phosphorene sheets. The identical slope of the lines enforces consistency in the extrapolation for calculations where the 
distance between adjacent sheets are varied. 

This difficulty has been noted in previous quantum Monte Carlo studies of layered materialPS where a combination 
of the KZK functional 810 was used for correlations in the plane and an extrapolation scheme was used for correlations 
in the direction perpendicular to the planes. Our approach is to extend this methodology in two regards. Firstly, we 
utilize extrapolation in the size of the supercell to determine correlations both in plane and out of plane and secondly, 
because such extrapolations can introduce a significant amount of noise in the extrapolated quantities, we correlate 
the parameters of the extrapolations between systems with similar geometries. This is best illustrated looking at the 
extrapolation of the energy in the direction perpendicular to the phosphorene planes. 

This correlation may be expected to decrease with distance between the planes and also be a function of the number 
of layers simulated. To capture this, we extrapolate based on the total distance between images in the perpendicular 
direction and also require that the slope of the extrapolation be a simple function of the distance between the planes. 
The results of this procedure are shown in fig |ST| where the energy per atom for different supercells with different 
spacings between the phosphorene planes and different numbers of layers in the supercell are shown as a function of 
the inverse of the distance between the top and bottom of the supercell. All of these calculations had a three by three 
tiling of the primitive cell of phosphorus in the direction parallel to the planes and both the spacing between the 
layers and the number of copies of the supercell were varied. A similar procedure was used to extrapolate to infinite 
size supercells in the other direction. 
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II. NUMERICAL CONSIDERATIONS 


To control computational cost, we verify that our choice of DMC time step satisfies a balance between efficient 
sampling and having a large time step bias. In doing so, for each system we choose fixed moderately sized supercells 
and perform short DMC runs at different time steps. A best fit line is constructed for the energy as a function of 
time step and we choose to use the largest time step that is within 1 mHa/atom of the zero time step extrapolated 
energy. We have found that a time step of 0.0075 a.u. is adequate to achieve this level of accuracy in all cases. 

Further, to control the memory required by the wave function we assess the effect of enlarging the grid spacing 
associated with the B-spline representation of the Kohn-Sham orbitals relative to the equivalent real space grid used in 
the plane wave pseudopotential calculation in which they are generated. The convergence of the total energy, kinetic 
energy, and variance in the local energy are assessed in determining the appropriate grid spacing. Noting that the 
variance converges most slowly in the grid spacing, an enlargement factor of 4/3 preserves accuracy while reducing 
the memory used up in representing the wave function. 


III. PSEUDOPOTENTIAL GENERATION AND TESTING 

Two different pseudopotentials were generated for this study, one treating 5 valence electrons and the other treat¬ 
ing 13. Both pseudopotentials were generated using the opium pseudopotential generation codeP^ The benchmark 
quantities under consideration are the equilibrium properties of a phosphorus dimer and the ionization potential and 
electron affinity of an isolated phosphorus atom. 

In Table [I] some atomic benchmarks are given. While the accuracy relative to experiment is comparable for the 
ionization potential for both pseudopotentials, the 5 electron pseudopotential performs worse for the electron affinity. 
Given the relatively primitive trial wavefunction used, these results should not be taken as conclusive, but errors due 
to the nodal surface are generally less than a few tenths of an eV. 



5 electron 

13 electron 

Experiment 

Ionization Potential 

10.7112 ± 0.00084 

10.6832 ± 0.0598 

10.48669 

Electron Affinity 

0.6405 ± 0.0084 

0.7483 ± 0.0626 

0.746609 


TABLE I. Ionization potential and electron affinity of an isolated phosphorus atom computed using two different pseudopoten¬ 
tials. 

We also calculate the binding energy, atomization energy and vibrational frequency for the phosphorus dimer by 
calculating the energy as a function of bond length. The results of these calculations for each pseudopotential and a 
fit to a morse potential are shown in Figure |S2| and the resulting observables are summarized in Table \H\ 



FIG. S2. Phosphorus dimer energy vs bond length curves for 5 and 13 electron pseudopotentials. Experimental values for the 
ionization potential and electron affinity are taken from References |S2'1 and IS3l respectively. 
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5 electron 

13 electron 

Experiment 

bond length (A) 

1.8618 ± 0.009 

1.8824 ± 0.0018 

1.89340 ± 0.00044 

vibrational frequency (cm -1 ) 

827.3 ± 9.4 

851.7 ± 31.8 

780.77 

atomization energy (eV) 

4.900 ± 0.001 

4.894 ± 0.040 

5.03 ± 0.02 


TABLE II. Properties of phosphorus dimer computed using two different pseudopotentials. The experimental bond length, 
vibrational frequency and atomization energy are taken from references IS4l [Sol and EU respectively. 


While both the 5 electron and 13 electron pseudopotentials have similar accuracy relative to experiment, calculations 
using the potential with 13 electrons in the valence are nearly 100 times more expensive than those using the 5 electron 
potential. Therefore, its success in these metrics leads us to use the 5 electron pseudopotential for all remaining 
calculations on bulk phosphorus and phosphorene. 
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